(scl 4)
(seed (in "/dev/urandom" (rd 8)))

(de fact (N)
   (if (=0 N) 1 (apply * (range 1 N))) )
            
(de analytical (N)
   (sum
      '((I)
         (/ 
            (* (fact N) 1.0) 
            (** N I) 
            (fact (- N I)) ) )
      (range 1 N) ) )

(de testing (N)
   (let (C 0  N (dec N)  X 0  B 0  I 1000000)
      (do I
         (zero B)
         (one X)
         (while (=0 (& B X))
            (inc 'C)
            (setq
               B (| B X)
               X (** 2 (rand 0 N)) ) ) )
      (*/ C 1.0 I) ) )

(let F (2 8 8 6)
   (tab F "N" "Avg" "Exp" "Diff")
   (for I 20
      (let (A (testing I)  B (analytical I))
         (tab F
            I
            (round A 4)
            (round B 4) 
            (round 
               (*
                  (abs (- (*/ A 1.0 B) 1.0)) 
                  100 ) 
               2 ) ) ) ) ) 

(bye)
   
